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INTRODUCTION 


The  successful  development  of  an  autonomous  vision  system  for  mobile 
vehicles  would  be  of  considerable  value  and  importance  to  defense  and  related 
fields.  Numerous  reports  and  studies  currently  recommend  artificial 
intelligence/robotics  applications  which  require  autonomous  vehicles. 

Essential  to  these  robotic  vehicles  Is  an  adequate  and  efficient  computer 
vision  system.  A  potentially  more  successful  approach,  other  than  TV  pictures 
and  photographs,  would  be  to  develop  a  three-dimensional  system  employing  a 
laser  rangefinder. 

A  range  matrix  describing  a  certain  scanned  area  of  the  terrain  in  front 
of  the  mobile  robot  can  be  used  to  estimate  the  slopes  of  the  terrain.  The 
in-path  and  cross-path  slopes  of  the  terrain  are  evaluated  by  a  slope 
estimation  scheme.  These  slope  informations  along  the  passible  corridors  are 
utilized  to  determine  a  safer  and  more  accurate  path  for  the  mobile  robot 
vehicle  to  travel. 

The  mobile  robot  vehicle  is  equipped  with  data  acquisition  and  decision 
making  devices  for  its  autonomous  navigation  over  rough  terrain.  A  laser 
rangefinder  can  be  operated  by  emitting  laser  pulses  and  measuring  the  time  of 
flight  of  a  pulse  between  the  instant  it  is  transmitted  and  the  instant  the 
reflected  pulse  is  received.  This  time  of  flight  is  related  to  the  distance 
between  the  transmitter  and  the  point  on  the  terrain  from  which  the  pulse  is 
reflected.  The  terrain  is  scanned  by  changing  the  azimuth  and  elevation 
angles  of  the  laser  beam  in  a  discrete  fashion.  The  measurements  are  then 
available  in  the  form  of  an  NXM  'range-matrix'. 

References  are  listed  at  the  end  of  this  report. 


The  slope  estimation  problem  dealt  with  here  is  one  of  obtaining  smoothed 
estimates  of  function  values  and  particularly  their  derivatives  from  a  finite 
set  of  inaccurate  measurements  in  two  dimensions.  In  one  approach  we  can 
identify  the  dynamic  equations  of  the  underlying  system,  or  estimate  the 
distributions  for  the  quantities  of  interest  and  then  apply  optimal  estimation 
algorithms.  In  some  engineering  problems  the  stochastic  system  may  not  be 
identified  easily  and  in  these  situations,  spline  smoothing  has  proved  to  be  a 
useful  alternative. 

In  this  report,  we  obtain  the  smoothed  estimates  of  the  slopes  by 
utilizing  a  two-dimensional  smoothing  algorithm.  For  the  problem  of  smoothing 
a  finite  set  of  noise  corrupted  data  of  an  unknown  function,  it  is  proposed  to 
obtain  the  smoothed  estimate  by  fitting  a  two-dimensional  approximating 
function  to  the  data  set,  for  a  set  of  measurements  corrupted  by  a  white  noise 
process. 

HISTORICAL  REVIEW 

By  noting  the  fact  that  original  signals  such  as  visual  scenes  are  in 
analog  form,  techniques  were  developed  which  reconstruct  analog  signals  from 
discrete  data  by  utilizing  interpolation  or  approximating  functions. 

Frequency  domain  interpretation  of  the  interpolation  process  was  reported  in 
Reference  1.  Also,  B-spline  interpolates  (refs  2-4)  were  used  (ref  5)  in 
restoring  a  continuous  signal  from  a  set  of  digitized  data.  For  one¬ 
dimensional  noise  corrupted  data  generated  by  unknown  systems,  Reinsch  (ref  6) 
utilized  natural  cubic  splines  (refs  2-4)  along  with  least  squares  constraints 
to  solve  the  problem  of  curve  plotting.  Hou  and  Andrews  (ref  7)  constructed 
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continuous-discrete  image  and  utilized  spline  basis  functions  along  with  the 
least  squares  constraints  for  image  restoration.  Because  of  their  non¬ 
recursiveness,  the  algorithms  in  References  6  and  7  are  involved  with  complex 
computations  and  cannot  be  implemented  on-line.  Recently,  by  using  a 
reproducing  kernel  Hilbert  space  approach,  Weinert  et  al  (refs  8,9),  developed 
a  structural  correspondence  between  spline  interpolation  and  linear  least 
squares  smoothing  of  a  particular  random  process. 

In  recent  years,  two-dimensional  recursive  filters  have  drawn  much 
attention  because  of  the  need  for  processing  images  or  other  two-dimensional 
information.  Previous  efforts  (refs  10-12)  to  achieve  a  truly  recursive 
two-dimensional  filter  were  of  only  limited  success  because  of  the  difficulty 
in  establishing  a  suitable  two-dimensional  recursive  model  as  well  as  the  high 
dimension  of  the  resulting  matrix  and  state  vector.  Recently,  by  using  a  two- 
dimensional  recursive  model  obtained  from  a  two-dimensional  spectral 
factorization  technique  (ref  13),  Woods  and  Radewan  (ref  14)  developed  a  two- 
dimensional  Kalman  vector  processor  and  a  two-dimensional  Kalman  scalar 
processor. 

The  above  mentioned  time-domain  design  techniques  assumed  or  identified  a 
two-dimensional  stationary  discrete  system  model  at  the  beginning  of  their 
problem  formulation.  On  the  other  hand,  Reinsch  (ref  6)  interpreted  a  one- 
dimensional  data  smoothing  problem  as  an  optimal  curve-fitting  problem  arising 
in  approximation  theory,  and  proposed  a  nonrecursive  smoothing  algorithm  using 
smoothing  splines.  For  a  two-dimensional  image  restoration  problem,  Hou  and 
Andrews  (ref  7)  followed  the  approach  taken  by  Reinsch  (ref  6),  and  extended 
it  to  a  two-dimensional  problem  in  a  nonrecursive  manner.  In  this  report,  we 
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develop  a  two-dimensional  recursive  smoothing  algorithm.  Compared  to  its 
nonrecursive  counterpart,  this  recursive  algorithm  will  require  less 
computational  complexity  and  memory  space.  Especially,  the  amount  of 
computation  needed  at  each  iteration  is  independent  of  the  size  of  the 
two-dimensional  data. 

PROBLEM  FORMULATION  FOR  ONE-DIMENSIONAL  APPROXIMATION 

From  the  viewpoint  of  approximation  theory,  when  a  set  of  discrete 
observation  data  is  noise  free,  spline  interpolation  provides  a  means  of 
optimally  reconstructing  an  unknown  original  signal.  When  the  observation 
data  are  corrupted  by  noises,  and  if  the  form  of  the  original  continuous 
signal  is  known,  then  we  can  use  least  squares  estimation  techniques  to 
approximate  the  original  signal.  We  are  now  dealing  with  a  problem  in  which 
an  unknown  signal  is  approximated  by  smoothing  splines  from  a  set  of  noise 
corrupted  observation  data.  Specifically,  an  unknown  signal  f(  5)  is 
approximated  by  a  polynomial  spline  s(  £)  which  minimizes  the  objective 
function: 

N  N 

J*  -  I  [*(5n)  -  an]TRn”l[8(5n)  “  «%,]  +  i  l  Pn L  [sk(5)]2dO  (1) 

n-1  «-2  V-l 

where 

mn  is  observation  data; 

®n  *  f(5n)  +  vn,  for  n  -  1,  2,...,N; 

vn  is  a  white  observation  noise  process  with  error  covariance 
A 

*n  "  E<vn*vnT>; 

Pn  >  0  is  a  smoothing  parameter;  and 
sk  is  the  kth  derivative  of  s(5). 
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At  this  point,  it  is  worthwhile  to  note  the  physical  role  of  the  smoothing 
parameter  pn  as  follows:  (1)  when  pn  becomes  very  small,  pn  0+,  the 
resultant  approximating  function  will  pass  through  each  data  point  and  become 
an  interpolation  function;  (2)  when  pn  assumes  a  very  large  value,  pn  -*■  °°, 
minimization  of  the  objective  function  in  Eq.  (1)  corresponds  to  fitting  a 
straight  line  to  a  data  set  using  least  squares  criterion.  Thus,  it  can  be 
said  that  the  smoothing  parameter  controls  resolution  in  a  tradeoff  of  the 
smoothness  of  the  restored  function. 

Choice  of  Approximating  Function 

As  has  been  noted,  it  is  desired  to  develop  a  recursive  algorithm  whose 
results  are  sufficiently  close  to  those  obtained  by  directly  minimizing  the 
global  problem  as  given  by  the  criterion  in  Eq.  (1).  Fundamental  problems 
encountered  in  developing  a  recursive  algorithm  which  generates  approximating 
functions  are: 

1.  feasibility  of  recursive  structures, 

2.  feasibility  of  numerical  calculations. 

Regarding  the  first  problem,  it  has  been  noted  from  References  1  through  5 
that  some  of  the  approximating  functions  such  as  polynomial  splines  and 
piecewise  Hermite  polynomials  have  finite  support.  That  is,  a  resultant 
approximating  function  for  one  section  is  mostly  affected  by  its  neighboring 
data  points.  Thus,  a  recursive  structure  with  one  or  more  sample  delays  would 
result  in  sufficiently  close  results  to  nonrecursive  ones. 

For  the  second  point,  out  of  a  certain  set  of  functionals,  an  optimal 
solution  to  Eq.  (1)  is  an  L-spline  (refs  2-4).  L-spline  is  a  piecewise 
polynomial  of  degree  2k-l,  and  has  2k-2  continuous  derivatives  in  the  region 


5 


[Cl,  ?n3 •  Here,  we  propose  to  restrict  our  approximating  functions  to 
piecewise  Hermite  polynomials  (ref  3)  of  degree  2k-l ,  which  have  k-1 
continuous  derivatives  in  the  region  ICi,Cn]»  Advantages  in  using  piecewise 
Hermite  polynomials  are  as  follows.  Define 


A  .  .  I 

-  [s(C),  s'(C),...,sk_1(0]T|  ,  i  =  1,...,N 

I 

Then  a  piecewise  Hermite  polynomial  s( Q  is  completely  determined  by  xif  i  = 
l,2,...,N.  For  the  purpose  of  clarity  in  discussion,  only  the  case  of  k=2  is 
treated  in  the  following.  A  piecewise  cubic  Hermite  polynomial  is  represented 
as: 

<  < 

sl,2(  5)  for  =  £  -  C2 

•  • 

s(S)  -  .  .  (2) 

<  < 

SN-1 ,N'  for  %-l  *  5  «  % 

where 


sk-l,k(5)  “  [<t>k,l(5)4k,i(0<l>k,o(0  'Kt,o(5)][xkT,  x^1]1 
xk  “  [s(5k)  s'(Ck)]T,  k  -  1 ,  • . .  ,N 
*k,l(5)  *  (?-Sk-l>2I(5k-Sk-i)  +  ZCCk-Ol/CCk-Cfc-!)3 

«M(C)  -  (s-Sk-i^ce-^/cCk-Ck-!)2 
4>k,0(c)  “  (5k-?)2t(5k-^k>l)  +  2(C-Ck-1)]/(Ck-Ck-1)3 
*k,0<«)  *  C^^-xKCk-Q^CCk-Ck-i)2 


Smooth  Integrals  as  Quadratics  at  Node  Points 


Thus,  it  becomes  natural  that  the  smoothing  integral  in  Eq.  (1)  is 


(2a) 
(2b) 
(2c) 
(2d) 
(2e) 
(  2  f ) 


expressed  in  terms  of  xj/s,  i  -  1,2,...,N.  With  some  manipulations  in 
algebra,  the  smoothing  integral  for  k  -  2  in  Eq.  (1)  is  represented  in  a 
quadratic  form  as  derived  in  Appendix  A. 


where 


/n  !  I  s"(  C)  I  I  2  dC  =  (Xn  -  A*xn.1)TB-1(xn-A*xn_1) 
5n-l 

l*n-l  T  CU  C12  Ixn-i 

*  I  •  *1 

j*n  C2i  C22  jxn 


[s(Q,  s'(Q] 


|1  A| 

A*  -  |  !  B 

10  1| 


(3) 


(4) 


(5a) 


A  =  5n  “  ^n-l.  for  n  =  2,...,N  (5b) 

A  A  A  A 

Cu  -  A*tB_1A*  ,  C12  -  C21t  -  -A*T  B"1,  C22  -  B-1  (5c) 

A  nonrecursive  solution  for  the  minimization  of  the  objective  function  in 
Eq.  (1)  can  be  obtained  by  taking  the  gradient  of  J*  with  respect  to  [x^, 
x2,...,Xjj]T  and  setting  it  to  zero.  However,  this  approach  will  require 
solutions  of  a  set  of  2N  simultaneous  equations.  To  avoid  this  computational 
problem,  we  have  developed  a  recursive  algorithm  which  requires  inversion  of  2 
x  2  matrices  only. 

RECURSIVE  ALGORITHM 

Given  a  set  of  initial  values  for  the  mean  and  its  error  covariance 
P^,  where  =*  E{(x^-X|)(x^-x^)T },  by  using  Eqs.  (3)  and  (4),  the  objective 
function  in  Eq.  (1)  becomes 
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discrete  edges,  a  two-dimensional  smoothing  algorithm  is  utilized  to  estimate 
the  range  slopes;  (3)  estimated  range  slopes  are  transformed  into  terrain 
slopes. 

Estimated  Terrain  In-Path  Slope 

The  simulation  of  terrain  with  hills  and  valleys  is  given  in  Figure  2. 

The  estimated  terrain  in-path  slopes  (ref  18)  are  displayed  in  terms  of  a 
slope  map,  Figure  2.  Characters  A,...,G  represent  a  particular  range  of  the 
terrain  in-path  slopes  increasing  from  A  to  G,  at  the  corresponding  location. 

U  represents  undefined  slopes.  In  Figure  3,  we  note  circular  slope  regions  on 
the  faces  of  sinusoidal  hills  and  valleys.  Also,  along  a  radial  direction, 
the  estimated  slopes  are  changing  slowly  from  one  region  of  slopes  to  another. 
The  large  empty  spaces  are  due  to  the  hidden  regions  at  the  back,  of  boulders 
or  hills  where  laser  rays  could  not  reach.  The  undefined  gradient  represented 
by  'U'  occurs  when  the  recursive  algorithm  cannot  be  applied  due  to  sharp 
changes  in  ranges  between  adjacent  measurement  data.  The  estimated  in-path 
terrain  slope  maps  are  used  for  the  evaluation  of  the  terrain  In  front  of  the 
mobile  robot  vehicle. 

Terrain  Cross-Path  Slopes 

In  discussing  the  terrain  cross-path  slopes  (ref  19),  the  data  can  be 
conveniently  processed  to  generate  smoothed  in-path  and  cross-path  range 
slopes  recursively  in  a  spherical  coordinate  system  due  to  the  regularity  of 
the  elevation  and  azimuth  angles.  When  we  proceed  to  calculate  the  true 
terrain  slopes  on  the  base  plane,  the  regularity  of  the  data  points  is 
completely  destroyed.  For  a  fixed  elevation  angle  S,  the  horizontal 
projection  of  the  range  data  is  not  located  at  a  fixed  distance  from  the 
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measurement  data  tn  cylindrical  coordinates.  However,  there  is  a  major 
difficulty  in  this  approach.  Even  though  the  two  independent  variables  8^  and 
9 j  for  the  rangefinder  are  changing  with  constant  increments  A8  and  A0, 
respectively,  the  independent  variable  in  a  cylindrical  coordinate  changes 
irregularly.  The  recursive  smoothing  algorithm  in  the  previous  subsection 
requires  that  the  data  points  be  located  at  the  corners  of  rectangular  grids 
of  the  two  independent  variables.  Since  the  two  independent  variables  p^  and 
9j  in  a  cylindrical  coordinate  system  do  not  form  rectangular  grids,  the 
smoothing  algorithm  cannot  be  applied  directly.  By  noting  that  the 
positioning  angles  8^  and  9j  are  changing  in  regular  fashion,  it  is  proposed 
to  obtain  the  smoothed  estimates  of  the  range  slopes  dr/d  8  and  dr/d  9  defined 
in  spherical  coordinates.  Then,  these  estimates  are  transformed  to  the 
terrain  slopes.  In  applying  the  smoothing  algorithm  to  terrain  slope 
estimation,  one  point  to  be  mentioned  is  that  the  basic  philosophy  of  the 
smoothing  spline  approach  is  to  suppress  the  noise  elements  by  fitting  a 
smooth  approximating  function  to  a  noise  corrupted  data  set.  Therefore,  when 
the  function  to  be  approximated  has  sharp  changes  in  its  values  or 
derivatives,  the  smoothing  algorithm  will  produce  errors  in  the  results  by 
smoothing  out  these  actual  sharp  changes.  From  the  viewpoint  of  terrain  slope 
estimation,  such  changes  occur  at  the  edges  of  a  boulder,  a  crater,  or  a  ridge 
on  the  terrain.  Thus,  it  is  proposed  to  detect  these  edges  by  using  the  rapid 
estimation  scheme.  Then,  for  the  area  which  is  free  of  discrete  edges,  the 
two-dimensional  smoothing  algorithm  is  utilized  to  estimate  the  slopes.  The 
terrain  slopes  are  estimated  in  the  order:  (1)  discrete  edges  are  detected 
by  using  the  rapid  estimation  scheme;  (2)  for  the  area  which  is  free  of 
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in  the  fourth  column,  x*k,4|k,4  and  x*k,4lk+l,4»  an<*  80  on* 

The  approximation  method  described  above  is  one  of  the  simplest  ones*  We 
can  employ  more  elaborate  approximation  methods  at  the  cost  of  more 
complicated  computations.  By  now,  we  have  introduced  a  recursive  quarter- 

A 

plane  processor  which  computes  the  filtered  estimate  x*^  £|k  £  as  an 

A  A 

approximation  to  xk,Z|k,i*  From  the  definition  of  the  estimate  xk,il|p,q>  the 
value  of  ^  jj  which  minimizes  J  is  Xk  £|n  m*  Here,  we  note  that  x*k^  £.|  k,  X.  Has 

A 

its  support  in  the  region  R-(k,&)»  while  the  nonrecursive  solution  x*k,Jl|N,M 
has  its  support  in  the  region  R({j  j,j)*  Thus  if  we  desire  to  have  a  better 

A 

approximation  to  x^  £|jj  we  need  to  develop  a  smoothing  algorithm  which 

>  <  < 

computes  x^  £ | k+d  £+<j  where  d^,d2  ■  1,  k+d^  *  N  and  H+d2  ■  M. 

*  1  ’  2  . 

The  smoothed  estimate  x^  jj,|k+d  l+d  *s  defined  as  the  estimate  of  x^  £ 

*  1  2 

obtained  by  fitting  an  approximating  function  in  the  region  R(k+dq ,  * 

The  derivation  procedure  for  the  above  approximation  is  similar  to  that  of 
filtering  discussed  previously,  and  is  omitted  for  conciseness. 


FURTHER  NAVIGATION  PROBLEMS 
Terrain  Slopes  and  Range  Slopes 

With  reference  to  Figure  1,  terrain  in-path  and  cross-path  slopes  are 
defined  as  the  two  orthogonal  slopes  dz/dp  and  dz/pd6  in  a  cylindrical 
coordinate  system.  During  the  past  investigations,  the  terrain  slopes  were 
found  to  be  appropriate  measures  for  evaluating  a  terrain.  A  direct  approach 
for  estimating  the  terrain  slopes  would  be  to  fit  a  smoothing  spline  to  the 


By  definition,  the  filtered  estimate  x^  2|k,2  is  the  estimate  of  x^^ 
obtained  from  an  approximating  function  which  minimizes  the  objective  function 


in  Eq.  (27)  for  (p,q)  =*  (k,2).  For  minimization,  we  take  the  gradient  of 


J(k,2)  with  respect  to  xi  and  j£2»  and  set  it  to  zero 


where 


V(*l)  J(k,2)  "  0 
*2 


Xj  ”  t*l,jT.  xzj1 


.Xk,jT]' 


We  can  obtain  a  final  recursive  estimation  equation  as: 


xk,l |k,2 


Px.r1  xk,  1 


xk-l ,  1 1  k-1 , 2 


xk,2|k,2 


aTRk,2"1  mk,2 


xk-l, l|k-l,2 


For  notation  used  in  the  above  equation,  see  Reference  17.  With  reference  to 


the  final  estimation  equation  above,  it  is  noted  that  for  the  filtered 


estimate  x^  2|k  2  the  scheme  uses  the  previous  estimates  xk  3,  x^_^  1 1  k.— 1  2* 

A 

and  Xk_i  2 1 k-1  2>  and  the  measurement  m^  2*  Here,  it  should  be  emphasized 

A  A 

that  the  scheme  uses  the  smoothed  estimate  Xfc_i  1 | k-1  2  instead  of  x^-i  ^ . 

A 

Thus,  when  the  filtered  estimate  xj^  2 ) k  2  *-s  computed,  we  need  to  update  the 

A  A 

estimate  x^i  to  %  2  ^or  u8e  *-n  next  iteration.  Now,  by  using  the 


recursive  estimation  equation  in  Eq.  (28)  and  the  pseudo  error  covariance 


equation  in  Reference  17,  we  can  compute  x^  2 | k  2»  k  *  2,...,N  recursively. 

A 

The  resultant  recursive  filtering  equation  for  x^^jk  3  becomes  similar 

A 

to  the  one  in  Eq.  (28).  Also,  the  smoothing  equation  for  x*i  3/i+l  3  becomes 

A 

similar  to  the  one  for  x*^  2/i+l  2*  After  all  the  estimates  in  the  third 

A  A 

column,  x*i  3|k>3  and  x*^  3 1 3  for  k  ■  1,  2,...,N  are  obtained,  the 

A 

smoothed  estimates  x*^ 31^+1^,  k  *  1,...,N-1,  will  be  used  for  the  estimates 


r-  **-  f.f.  .. 


A  QUARTER-PLANE  PROCESSOR 

A 

The  estimate  x^  £|p  q  is  defined  as  the  estimate  of  obtained  by 

fitting  as  approximating  function  in  the  region  R(p,q)  where 

A  <  <  <  < 

R(p,q)  "  {(S,h)l  *1  “  5  •  h  »  and  i*  -  n  «  V  (26) 

>  > 

For  (p,q)  *  (k,A),  x^  £|k  £  becomes  a  filtered  estimate.  For  p  *  k  and  q  =  £, 

A 

except  for  p  »  k  and  q  —  £,  xk,A|p,q  becomes  a  smoothed  estimate  of  xj^  £•  In 

A 

our  formulation,  xjt>£|Pjq  would  be  obtained  by  minimizing  the  objective 
function  in  the  region  R(p,q): 

J(P»q)  =  I  I  [(Hxi  j-mi  j)^”1  (Kxi  j-mi  j)] 

1-1  i-1  i,J  J 


+  [  l  l  <xi,jT»  xi+l, jT»  xi , j+lT »  xi+l , j+lT> 
j-1  i=l 

•  C(xltjT,  x1+ijT,  Xi,j+1T»  xi+l,j+lT)Tl 
where  H  -  (1,0, 0,0). 


(27) 


In  a  quarter-plane  processor,  the  filtered  estimate  xk>£|kj£  is  obtained 
by  using  the  previous  estimates  of  x^-i  ^  £_^ ,  xjt-i  £,  and  x^^-i,  and  the 
measurement  In  the  next  iteration,  the  filtered  estimate  x^+i  >  £ |  fc+j.  is 

obtained  by  using  the  previous  estimates  of  Xfc  £-1,  £,  and  x^+i  >  £-.^  and  the 

measurement  m^+i  £.  After  estimating  all  the  states  in  the  ith  column  the 

A 

recursive  processor  moves  to  the  next  column,  and  estimates  ^  JJ.-+-1 1  k. ,  £-t-l » 
k  -  2,...,N,  and  so  on.  First,  we  will  discuss  a  filtering  procedure  for 


xk,2|k,2*  k  *  2,  3,...,N.  Then,  this  procedure  is  extended  to  the  filtered 
estimates  x^  £1^  £  for  l  *  2,  3,...,M. 


so-called  "saddle  point"  and  every  surface  element  of  such  a  membrane  is  "pure 
twist."  An  appropriate  smoothness  measure  would  be  changed  to: 


z(5,h)  -  [--2  s(£,n)]2  +  [“2  sU.n)]2  (23) 

3.  In  Reference  7,  Hou  and  Andrews  suggested  using  I  I V4s( £, n) I  I  2  as  a 
smoothness  measure  for  a  surface.  The  physical  interpretation  of  the  quantity 
V4  s(£,n)  is  found  in  a  plate  bending  theory  (ref  16);  an  unloaded  plate  can 
bend  only  in  a  biharmonic  function  w  where 

V4w  -  0  (24) 

A  bicubic  Hermite  polynomial  which  minimizes  the  objective  function  J  is  given 
in  Appendix  B. 

Smoothing  Integral 

Now,  it  is  needed  to  determine  the  function  s(£,n)  which  minimizes  the 

objective  function  J  in  Eq.  (16).  It  is  noted  that  the  smoothing  integral  in 

its  present  form  gives  difficulties  in  finding  an  explicit  solution.  By 

evaluating  the  Integrals  of  the  derivatives  of  basis  functions  and  applying 

some  algebraic  manipulations,  these  smoothing  integrals  are  converted  to 

quadratic  forms  as  follows: 

hM  eN  M-l  N-l  n-j+i  Ci+i 

p  /  /  [z(C,h)]d£dn  -Ilf  f  [z(£,n)]d£dn 

hi  e 1  j-l  hj  £i 

“  X  \  (Xi’JT’  Xi+1*JT’  x2+l,j+lT)* 

j“l  i-1 

c*(xi,jT,  xi+1>jT,  xi>j+1T,  xi>;j+1T)T 
where  C  is  a  16  by  16  matrix. 


(25) 


for  i  »  1,  2,  ...,N  and  j  -  1,  Then,  a  piecewise  bicubic  Hermite 

polynomial  is  completely  defined  by  x^j,  i  -  1,  2,...,N  and  j  ■  1,  2 . . 


as  follows: 
and 

where 


<  < 

s(5,n)  -  stjCC.n)  ,  for  -  5  -  5i+i 

<  < 


1  1 

*i,j(5.n)  -  I  l 

£-0  m-0 


n  -  rij+1 

<t>£(  S)  *  <t\n(  n) 
♦t(  5)  •  4tn(  T1> 
♦t(  5)  *  n) ! 
♦t(  5)  *  4(  n) 


I  *  xi+A,j+m 


(18) 

(19) 


(20) 


Choice  of  the  Smoothness  Measure  z(5,n) 

Here,  we  present  three  examples  of  the  smoothness  measures  and  compare 
their  physical  implications. 

1.  Gaussian  curvature:  In  Reference  15,  the  mean  curvature  of  a 
surface  at  (£,n)  is  defined  as: 

(0.5)V2s(  5,n)  (21) 


Noting  the  Euler's  theorem  (ref  16)  that  the  sum  of  two  curvatures  in 
perpendicular  directions  at  a  point  is  constant,  the  square  of  V2s(£,n)  in  gq. 
(21)  would  be  a  reasonable  measure  for  the  smoothness  of  a  surface: 


z(C,n)  -  (  ---j  s(C,n)  +  ---5  s(5,n)] 
35  3n 


(22) 


2.  A  variation  from  the  Gaussian  curvature:  With  reference  to  Eq.  (22), 


an  interesting  case  occurs  when  the  two  principal  curvatures  are  equal  and  of 
opposite  sign.  The  mean  curvature  in  this  case  is  zero.  This  is  the 


PROBLEM  FORMULATION  FOR  TWO-DIMENSIONAL  APPROXIMATION 

When  the  observation  data  are  noise  corrupted  and  the  underlying  system 
is  unknown,  it  is  proposed  to  approximate  the  original  signal  by  spline 
functions  which  minimize  a  certain  objective  function.  Thus,  from  a  set  of 
discrete  measurements  corrupted  by  white  noise  process  vj^j 

+  vi,j  »  i  2,...,N,  j  ■  1,  2 . N  (15) 

The  original  two-dimensional  signal  f(C,n)  defined  in  the  region  of  (5,n)  is 
approximated  by  a  spline  function  s(5,n)  which  minimizes  the  following 
objective  function: 

N  N 

J  =  ^  [s(Ci,nj)  -  j ] Tr± ,  j"  X[s(  5i,  hj)  -  m^j] 

.’TM  .% 

+  P  [/  J  z(5,n)d£dn]  (16) 

nl  Ci 

where  p  >  0  is  the  smoothing  parameter;  Rj^j  is  the  observation  error 
covariance;  and  z(5,n)  is  a  certain  smoothness  measure  of  s(£,n)  at  (5,n). 
Choice  of  an  Approximating  Function 

In  this  report,  we  are  interested  in  obtaining  smoothed  estimates  of 
function  values  and  the  first  derivatives  in  both  5  and  n  directions.  Here, 
we  propose  to  restrict  our  approximating  functions  to  piecewise  bicubic 
Hermite  polynomials  which  have  continuous  first  derivatives  in  both  £  and  n 
directions. 

Define 

A  3s  3s 

xi,j  “  [s(C,n)  ,  “  (C,n)  ,  --  (5,n), 

,J  35  3n 

32s  J 

----  (€,n)jT|(5>n)  -  (q,^)  (17) 

3?3h  |  J 
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mi  “  f(5i)  +  Vi  » 


i  -  1, . . . ,100 


(12) 


where  is  white  Gaussian  measurement  noise, 

Rt  =■  ECviVi7)  -  0.000025,  and  =  £n  -  5n-i  *  2ir/100  -  0.062832  (13) 

Function  values  and  the  first  derivatives  at  discrete  nodes  are  estimated 
from  the  measurements  m^,  i  ■  1,...,100  by  the  three  schemes  below: 

1.  Difference  quotients  method. 

2.  Recursive  smoothing  algorithm  with  £  ■  1:  Eq.  (10). 

3.  Nonr.ecursive  smoothing  by  cubic  splines  as  described  in  Reference  6. 
Table  I  shows  the  mean-square  errors  from  the  three  schemes  above  where 


1  100 

e0  “  TTT  l  (f(5i)  >  Xi|i+1(l))2, 

100  1-1 


(14) 


l  100 

el  -  7“  l  (f*(«l)  -  M|1+1(2))2 
100  iml 


TABLE  I.  MEAN-SQUARE  ERRORS 


e0 

E1 


Difference 

Quotients 


2. 5x10" 5 
4. 30x10“ 1 


Recursive  Smoothing 
by  using  Eq.  (10) 


0.8x10“ 5 
0.12x10" 1 


Nonrecursive  Smoothing 
by  Cubic  Splines,  Ref.  6 


0.57x10“ 5 
0.1007x10“  1 


1 


From  Table  I,  it  is  noted  that  both  smoothing  algorithms  are  successful 
in  reducing  the  error  in  the  estimated  states.  The  error  in  the  first 
derivative  is  decreased  by  more  than  10  db.  Moreover,  Table  I  shows  that  the 
performance  of  the  two  smoothing  schemes  are  comparable.  However,  it  should 
be  emphasized  that  the  recursive  algorithm  developed  in  this  report  is  much 
simpler  than  the  nonrecursive  spline  smoothing. 
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following.  The  smooched  estimate  of  obtained  by  fitting  cubic  splines  to 

A 

n^,  n  ■  2,  3,...,  i+i  with  the  initial  values  and  E^  ■  pj  is  the  same  as 
the  smoothed  estimate  of  x^  obtained  by  fitting  splines  to  n^,  n  ■  i+l,...,l+A 

a 

with  the  filtered  estimate  xjj^  and  E^. 

A 

A  recursive  procedure  to  obtain  the  smoothed  estimate  x^|^+^  is 
summarized  as  follows: 

A 

Part  1.  Obtain  the  filtered  estimates  x^|^,  i  ■  2,...,N  by  using  Eqs. 
(7d) ,  (9),  and  (9a). 

A 

Part  2.  Obtain  smoothed  estimates  X£|i+i,  for  i  *  1,...,N-1  by  using 
Eqs.  (10),  (10a),  and  (10b). 

A 

As  was  mentioned  earlier  in  this  section,  the  smoothed  estimate  xjJi+jj,  is 

A  A 

an  approximation  to  the  nonrecursive  solution  Xj*  *  X^|N.  Thus,  as  the  number 
of  delays,  A,  increases,  we  will  get  a  better  approximation  to  x^*.  For  A  * 

2,  3,...,  only  the  smoothing  part  is  modified  by  solving  the  simultaneous 
equations  in  Eq.  (7.)  with  k  -  i+A.  In  fact,  the  smoothing  algorithm  developed 
by  far  is  a  fixed-lag  smoothing  algorithm,  which  is  suitable  for  an  on-line 
implementation.  If  the  situation  does  not  require  an  on-line  implementation, 
we  can  also  derive  a  fixed-interval  smoothing  algorithm  with  observation  set 
M  -  {m^ , . . . ,mfj) .  In  this  case,  the  resultant  smoothed  estimates  become 
exactly  the  same  as  the  nonrecurslve  estimates. 

SIMULATION  RESULTS 

For  a  continuous  signal 

f(C)  -  sin(f-)  (11) 

Measurements  are  obtained  at  discrete  points: 
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With  reference  to  Eqs.  (7)  and  (8),  each  iteration  of  the  recursive 


filtering  algorithm  can  be  interpreted  as  fitting  a  cubic  polynomial  to  the 

a 

previous  estimate  *1-1 | j-i  and  the  present  measurement  m^  in  the  region 


C^i-l.Si] • 

Smoothing 

A 

A  smoothed  estimate  x^|^+£  is  defined  as  the  estimate  of  x^  obtained  by 
solving  the  minimization  problem  in  Eq.  (6)  with  N  *  i+i.  In  fact,  this  can 
be  interpreted  as  fitting  a  polynomial  spline  to  the  first  i+l  data  and 
obtaining  the  function  value  and  its  derivative  from  the  approximating 
function  at  the  node  i.  In  our  formulation,  this  corresponds  to  solving  the 
simultaneous  equations  of  the  same  form  as  Eq.  (7)  with  j  ■  1,  k  *  i+4,  and 
with  the  quantities  and  d^  defined  in  Eqs.  (7a)  and  (7b).  By  using  the 
same  reduction  method  as  before,  the  result  would  be  the  same  form  as  itself 
(Eq.  (7)  with  j  -  1  and  k  *  i+£). 

For  the  case  of  a  one-sample  delay,  x^^  is  obtained  by  solving  the 
simultaneous  equations  in  Eq.  (7)  with  j  ■  1,  k  ■  i+l,  which,  in  turn,  yields 


xi  I  i+l  3  viEi~1*i|i  +  Kimi+1  (1°) 

where 

Vt  -  [-PiC12(PiC22  +  HTRi+r^PiCzi  +  Gi]-1  (10a) 

Ki  *  "viPic12(  Pic22  +  HTR1+r1H)-lHTR1+1“1m1+1  (10b) 

and  E£  is  defined  as  before. 

Equation  (10)  is  the  desired  smoothing  algorithm,  in  which  the  smoothed 

A 

estimate  xjj^+i  is  obtained  by  updating  the  filtered  estimate  x^|£  with  the 
measurement  .  The  smoothing  procedure  described  above  implies  the 


for  xjl-  Thus,  by  eliminating  the  first  four  equations,  i.e.,  the  variable  x^, 

A 

Eq.  (7)  is  reduced  to  the  same  form  as  itself  with  j  ■  2,  k  ■  i,  and  x2|2  and 
E2  are  calculated  by  Eqs.  (7c)  and  (7d).  Note  that  the  quantity  Ej  defined  by 

A  A 

Eq.  (7d)  is  called  pseudo  error  covariance,  because  -  E  {xjc“,Xk|k)(xk“xk|k)T} 

cannot  be  computed  in  a  recursive  manner  directly. 

By  applying  this  reduction  method  repeatedly,  the  original  equations  in 
Eq.  (7)  are  reduced  to  the  same  form  as  themselves  with  j  ■  i-1,  k  ■  1  and 
every  x j | j  and  Ej  is  computed  by  Eqs.  (7c)  and  (7d)  recursively.  Solving  Eq. 

(7)  with  j  »  i*l ,  k  ■  i  in  terms  of  x^,  we  obtain  a  recursive  estimate 
algorithm  as 

A 

xi|i  “  Ei t HrR-i”’ ~  Pi-lc21Gi-l“ldi-ll  (8) 

Equation  (8)  is  rearranged  as 

A  A 

*i|i  “  EiHTRi_lmi  +  Fixi-l|i-l  (?) 

where 

Fi  “  “Pi-lEic2l( Pi-lGll+Ei-l” *)~ lEi-l~  *  (9a) 

and  Ejl's  are  computed  by  Eq.  (7d)  recursively.  Equation  (9)  above  is  the 

A 

desired  filtering  equation  which  computes  XjJj.  from  the  previous  estimate 

A 

xi-l|i-l  and  the  present  measurement  m^.  From  the  viewpoint  of  smoothing 
splines,  the  recursive  filtering  algorithm  can  be  interpreted  as  follows.  The 
estimate  of  Xj[  obtained  by  fitting  cubic  splines  to  the  measurement  data  u^, 

A 

n  *  2,  3,...,  with  the  initial  values  x^  and  E^  -  is  the  same  as  the 
estimate  of  x^  obtained  by  fitting  a  cubic  polynomial  to  the  measurement  m^ 

A 

with  the  initial  values  at  stage  i-1,  xi-i|i-i,  and  Ef_i.  In  fact,  the  above 
interpretation  comes  from  the  mathematical  derivations  in  Eq.  (7)  through  Eq. 
(9). 
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(7) 

where 

j  ■  1  and  k 

=*  i,  and  G-t  and  d<  are  defined  by 

A  . 

Gj  “  Pjcll  +  Ej" 

(7a) 

A  * 

dJ  “  Ej'  Xjlj 

(7b) 

*j!j  “  EjfHTRj-Vj  +  Pj-iC 2^-1“ 

-ii »  *i i i  ■ 

xi 

(7c) 

A 

Ej”1  *  p  j— lc22  +  -  (  Pj-lC2l)Gj_i" 

1(Pj-lc12>. 

Ei  -  P,. 

(7d) 

It  should  be  noted  that  the  matrix  on  the  left  side  of  Eq.  (7)  is  diagonally 
dominant  and  positive  definite.  Here,  we  are  interested  in  solving  Eq.  (7) 


9 


(6) 


JN  “  I  [(Hxn'ran)TRn‘1(Hxn-ffln)]  +  (xl“xl)Tpl“  H^-xi) 
n»2 


+  1  pn ( x n“ A*x rv- 1  > TB ” 

n*2 


where  H  =  (1  0) 


Let  the  solutions  to  the  above  optimization  problem  be  [x^*  ,X2* ,  •  •  •  ,xjj*]  . 

a 

If  Xp|q  *-s  defined  as  the  estimate  of  Xp  obtained  by  minimizing  Eq.  (6)  with  N 

A  A 

=■  q,  then  Xj*  can  be  written  as  xjjjj.  Here,  it  is  proposed  to  approximate  a 

A  A 

nonrecursive  solution  Xj*  -  Xi|N  by  Xi|i+£,  where  l  -  0  and  l  «  N.  As  has 
been  mentioned  before,  due  to  a  local  base  property  of  the  polynomial  splines, 

A 

the  smoothed  estimate  would  be  sufficiently  close  to  the  nonrecursive 

A 

solution  Xj*  for  l  ■  1,  2,  or  3. 

Filtering 

A 

From  the  definition,  a  filtered  estimate  xjj^  would  be  obtained  by 
minimizing  the  objective  function  in  Eq.  (6)  with  N  ■  i.  By  taking  the 
gradient  of  Jj  with  respect  to  [x^,  X2,...,xj[],  we  have  a  set  of  simultaneous 
equations  as  follows: 


u 


T 


Il'AJH  ULf! 


oe 


rover.  It  is  desired  to  calculate  the  cross-path  slope  at  point  ((5^,9j). 
However,  in  general,  Pi,j  *  Pi,j+1*  Thus,  the  true  cross-path  slope  is  not 
along  an  arc  connecting  points  Pj_  j  and  j+l • 

Our  algorithm  to  calculate  the  terrain  cross-path  slope  can  then  be 
summarized  as  follows: 

1.  Obtain  the  range  measurements. 

2.  Use  the  smoothing  algorithm  to  calculate  the  range  cross-path  slope. 

3.  Obtain  the  terrain  cross-path  slope  and  its  variance. 

Evaluation  of  Terrain  Variables 

As  mentioned  in  the  introduction,  in  the  evaluation  of  terrain  variables, 
we  only  use  the  slopes  at  the  spine  and  track  points.  The  reason  is  twofold. 
First,  only  a  minor  part  of  this  path  selection  scheme  needs  the  data  of 
elevation.  Second,  if  we  adopt  the  elevation  estimates  as  our  input  data,  we 
will  get  larger  errors  in  the  calculation  of  in-path  and  tilt  slope  terrain 
variables. 

The  terrain  variables  and  the  variances  together  with  the  corresponding 
explanations  are  listed  below  (ref  20). 

1.  In-Path  Terrain  Variables.  The  in-path  slope  terrain  variable  gives 
the  average  of  the  in-path  slopes  for  the  four  vehicle  wheels  at  each  section. 
This  variable  is  a  measure  of  the  risk  in  the  forward  direction. 

2.  Tilt  Slope.  Tilt  slope  terrain  variable  is  used  to  estimate 
excessive  cross-path  slopes  which  may  cause  the  vehicle  to  tip  over. 

3.  Obstruction  Height.  The  obstruction  height  is  calculated  for  six 
different  locations  at  each  discrete  section  of  the  terrain.  The  maximum 
value  is  then  chosen  as  representative  of  this  whole  section. 


In  deriving  the  formula  for  a  typical  obstruction  height,  we  use  a  third 
order  polynomial  to  approximate  the  terrain  elevation  in  each  location.  By 
differentiating  this  polynomial  with  respect  to  the  distance,  we  get  an 
expression  for  the  slope.  With  the  known  data  of  slopes  at  the  three  points 
substituted  into  this  expression,  we  can  determine  the  coefficients  of  the 
polynomial.  Using  this  polynomial,  we  can  then  find  the  obstruction  height  in 
this  direction. 

4.  Wheel  Deviation.  The  wheel  deviation  variable  describes  the  offset 
of  any  of  the  four  wheels  from  a  plane.  Wheels  on  any  three  track  points 
define  a  plane.  For  each  combination  of  three  wheels  touching  the  terrain, 
the  deviation  of  the  fourth  wheel  with  respect  to  this  plane  is  defined  as  the 
wheel  deviation. 

A  set  of  the  measurement  data  is  obtained  by  the  described  scanning 
scheme.  The  range  measurement  data  are  processed  by  the  gradient  estimation 
scheme  to  evaluate  ln-path  an!  cross-path  slopes  at  the  data  points.  Since 
the  slopes  are  estimated  in  the  spherical  coordinate  system,  it  is  necessary 
to  transform  the  range  slope  in  the  spherical  coordinate  system  to  terrain 
slope  in  the  cylindrical  coordinate  system.  The  in-p^th  and  cross-path  slopes 
and  their  covariances  at  the  spline  and  track  points  along  the  corridors  are 
evaluated  by  applying  a  two-dimensional  interpolation  scheme  over  the 
estimated  slopes  at  the  data  points.  The  terrain  variables  at  a  discrete 
section  along  each  corridor  are  computed  by  using  estimated  slopes  at  the 
spline  and  track  points.  Since  the  terrain  variable  estimates  have 
uncertainty,  the  present  method  increases  the  reliability  by  considering 
standard  deviation  as  well  as  mean  values. 


CONCLUSION 


By  taking  an  algebraic  approach,  a  recursive  smoothing  algorithm  was 
developed  as  an  approximation  to  nonrecursive  spline  smoothing.  Compared  to 
the  recursive  smoothing  algorithm  suggested  by  Weinert,  the  smoothing 
algorithm  in  this  report  is  simpler  in  that  the  scheme  is  in  a  discrete  form. 
Simulation  result  shows  that  the  performance  of  the  recursive  smoothing 
algorithm  is  comparable  to  that  of  its  nonrecursive  counterpart.  In  addition, 
the  computational  complexity  with  recursive  smoothing  algorithm  is  much  less 
than  its  nonrecursive  one.  Also,  recursive  smoothing  by  splines  can  be 
implemented  on-line.  By  taking  an  algebraic  approach,  a  two-dimensional 
recursive  smoothing  algorithm  was  developed  as  an  approximation  to  a 
nonrecursive  smoothing  spline  technique.  While  the  amount  of  computation 
required  for  a  nonrecursive  algorithm  increases  rapidly  with  the  size  of  the 
two-dimensional  data,  the  amount  of  computation  for  this  smoothing  algorithm 
Increases  only  linearly. 
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figure  2.  A  slops  aap  in  the  r-y  plane  for  the  in-path  terrain  slopes. 
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APPENDIX  A 


EVALUATION  OF  SMOOTHING  INTEGRALS 
From  Eqs.  (2a)  through  (2e)  of  the  text,  a  piecewise  cubic  Hermite 
polynomial  in  the  section  is  represented  as: 


si-l,i(£) 


I  ♦ii(C)  lT 
I  ’  I 

I  *1,1(5)  I 

j  *i,0(5)  j 

L*i,0(Oj 


I  s(  5i)  I 
I  s'(Ci)  i 
I  s(Ci-i)  I 
Ls’cq-iU 


(A-l) 


where  s(Si-i),  s'(Si-i),  s(Si),  and  s'(5i)  are  the  function  values  and  first 
derivatives  at  the  nodes  i-1,  and  i.  We  make  the  change  of  variables  such  as 

U  -  5  -  Si-i  (A-2) 

This  change  of  variables  does  not  affect  the  value  of  the  smoothing  integral 
and  results  in  a  simpler  computation.  The  smoothing  integral  in  the  interval 
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ii-1,1  -  /  lls"i_1.i<5)ll2  dC  -  /  Ms,,i-1,i(u)ll2du  (A-3) 

Z+l-l  J  0+  J 


where 

A  -  Ui  -  Ui-1  -  Si  -  Si-1 

Using  Eq.  (A-2)  and  Eq.  (2a)  of  the  text,  Eq.  (A-l)  becomes 

si— l ,i< -  [4>i,i(w)*i,i(y)  <t>i , o C u)  *i,o(  w)J  txiT  *i-iT]T  (A-4) 
Thus,  the  norm  square  of  the  second  derivative  is  written  as: 
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By  utilizing  Eq.  (A-5),  Eq.  (A-3)  becomes: 
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Thus,  smoothing  integral  is  obtained  as: 
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and  A*  as  below: 
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Equation  (A-7)  is  rewritten  in  the  following  form  as: 
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which  is  Eq.  (3)  in  the  text. 


APPENDIX  B 


In  this  Appendix,  we  show  that  a  piecewise  bicubic  Hermite  polynomial 
which  minimizes  the  objective  function  in  Eq .  (B-l)  becomes  a  bicubic  spline. 


J  -  JE  +  PJS 


where 


and 


M  N 


JE  3  ^  [s(?i,nj)  -  mlf  j]TRij~  L[s(  Ci,  Hj)  -  mij] 

dm  ?n  a4 

Js  ■  f  /.  (---2---2  s(C,n))2d?dn 


(B-l) 


(B-2) 


(B-3) 


Let  a  set  S  be  a  collection  of  all  piecewise  bicubic  Hermite  polynomials. 
Also,  we  define  a  set  U  as  a  collection  of  all  piecewise  bicubic  Hermite 
polynomials  which  satisfies  constraints  set  D  in  Eq.  (B-4). 

s(5i,nj)  *  c(i,j)  ,  i  *  1,2, ...,N  and  j  =>  1,2, ...,M 

=■  Cc(i,j)  ,  j  =  1,2, ...,M  and  i  -  1,N 

^(Sl.n-p/an  -  cn(i,j)  .  ,  1  -  1,2 . N  and  j  *  1,M 

923(5i,nj)/3C3n  -  C?>n(i,j)  ,  i  -  1,M  and  j  -  1,M  (B-4) 

Then  the  minimizing  problem  in  Eq.  (A-l)  is  rewritten  as: 


Min  J  ■  Min  [JE  +  pjg]  -  Min  [JE  +  p  Min  Js] 

s(c,n)es  s(C,n)cS  D  s(C,n)eU 


(B-5) 


In  the  paper  by  DeBoor  (ref  B-l),  it  is  noted  that  there  exists  a  unique 
bicubic  spline  g(C,n)  in  the  set  U.  Also,  by  using  a  standard  technique  to 
derive  the  minimum  norm  property  (ref  B-2)  of  a  bicubic  spline,  it  can  be 
shown  that: 


34g(C,n) 

J  /  ( — -—--)  2dCdn 

^1  5i  3?23n2 


(B-6) 
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Since  the  bicubic  spline  g(5,n)  is  unique,  we  have  the  following  Lemma: 

Lemma  1:  A  bicubic  Hermite  polynomial  s(5,n)eU  which  minimizes  the 
smoothing  integral  Jg,  becomes  a  bicubic  spline  g(£,n). 

With  reference  to  Eq.  (B-5)  and  Lemma  1,  we  conclude  that  a  piecewise  bicubic 
Hermite  polynomial  s(5,n),  which  minimizes  Eq.  (B-5),  becomes  a  cubic  spline. 


35 


APPENDIX  REFERENCES 


B-l.  DeBoor,  C.,  "Bibubic  Soline  Interpolation, "  J.  Math.  Phys.,  Vol.  41, 
1962,  pp.  212-218. 

B-2.  Ahlberg,  J.  H. ,  Nilson,  E.  N. ,  and  Walsh,  J.  L. ,  The  Theory  of  Splines 
and  Their  Application,  Academic  Press,  Inc.,  1967. 


TECHNICAL  REPORT  INTERNAL  DISTRIBUTION  LIST 

NO.  OF 
COPIES 


CHIEF,  DEVELOPMENT  ENGINEERING  BRANCH 

ATTN:  SMCAR-LCB-D  1 

-DA  1 

-DP  1 

-DR  1 

-DS  (SYSTEMS)  1 

-DS  (ICAS  GROUP)  1 

-DC  I 

CHIEF,  ENGINEERING  SUPPORT  BRANCH 

ATTN:  SMCAR-LCB-S  1 

-SE  I 

CHIEF,  RESEARCH  BRANCH 

ATTN:  SMCAR-LCB-R  2 

-R  (ELLEN  FOGARTY)  1 

-RA  I 

-RM  2 

-RP  1 

-RT  1 

TECHNICAL  LIBRARY  5 

ATTN:  SMCAR-LCB-TL 

TECHNICAL  PUBLICATIONS  &  EDITING  UNIT  2 

ATTN:  SMCAR-LCB-TL 

DIRECTOR,  OPERATIONS  DIRECTORATE  1 

DIRECTOR,  PROCUREMENT  DIRECTORATE  1 

DIRECTOR,  PRODUCT  ASSURANCE  DIRECTORATE  1 


NOTE:  PLEASE  NOTIFY  DIRECTOR,  BENET  WEAPONS  LABORATORY,  ATTN: 
OF  ANY  ADDRESS  CHANGES. 


SMCAR-LCB-TL, 


J7 


TECHNICAL  REPORT  EXTERNAL  DISTRIBUTION  LIST 


NO.  OF 
COPIES 

ASST  SEC  OF  THE  ARMY 

RESEARCH  &  DEVELOPMENT 

ATTN:  DEP  FOR  SCI  &  TECH  1 

THE  PENTAGON 

WASHINGTON,  D.C.  20315 

COMMANDER 

DEFENSE  TECHNICAL  INFO  CENTER 
ATTN:  DTIC-DDA  12 

CAMERON  STATION 
ALEXANDRIA,  VA  22314 

COMMANDER 

US  ARMY  MAT  DEV  &  READ  COMD 
ATTN:  DRCDE-SG  1 

5001  EISENHOWER  AVE 
ALEXANDRIA,  VA  22333 


COMMANDER 

ARMAMENT  RES  &  DEV  CTR 
US  ARMY  AMCCOM 

ATTN:  SMCAR-LC  1 

SMCAR-LCE  1 

SMCAR-LCM  (BLDG  321)  1 

SMCAR-LC S  1 

SMCAR-LCU  1 

SMCAR-LCW  1 

SMCAR-SCM-0  (PLASTICS  TECH  1 

EVAL  CTR, 

BLDG.  351N) 

SMCAR-TSS  (STINFO)  2 

DOVER,  NJ  07801 

DIRECTOR 

BALLISTICS  RESEARCH  LABORATORY  1 

ATTN:  AMXBR-TSB-S  (STINFO) 

ABERDEEN  PROVING  GROUND,  MD  21005 

MATERIEL  SYSTEMS  ANALYSIS  ACTV 

ATTN:  DRXSY-MP  1 

ABERDEEN  PROVING  GROUND,  MD  21005 


NO.  OF 
COPIES 

COMMANDER 
US  ARMY  AMCCOM 

ATTN:  SMCAR-ESP-L  1 

ROCK  ISLAND,  IL  61299 

COMMANDER 

ROCK  ISLAND  ARSENAL 

ATTN:  SMCRI-ENM  (MAT  SCI  DIV)  1 

ROCK  ISLAND,  IL  61299 

DIRECTOR 

US  ARMY  INDUSTRIAL  BASE  ENG  ACTV 
ATTN:  DRXIB-M  1 

ROCK  ISLAND,  IL  61299 

COMMANDER 

US  ARMY  TANK-AUTMV  R&D  COMD  1 

ATTN:  TECH  LIB  -  DRSTA-TSL 

WARREN,  MI  48090 

COMMANDER 

US  ARMY  TANK-AUTMV  COMD  l 

ATTN:  DRSTA-RC 

WARREN,  MI  48090 

COMMANDER 

US  MILITARY  ACADEMY 

ATTN:  CHMN,  MECH  ENGR  DEPT  1 

WEST  POINT,  NY  10996 

US  ARMY  MISSILE  COMD 
REDSTONE  SCIENTIFIC  INFO  CTR  2 

ATTN:  DOCUMENTS  SECT,  BLDG.  4484 

REDSTONE  ARSENAL,  AL  35898 

COMMANDER 

US  ARMY  FGN  SCIENCE  &  TECH  CTR 
ATTN:  DRXST-SD  1 

220  7TH  STREET,  N.E. 

CHARLOTTESVILLE,  VA  22901 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARMAMENT  RESEARCH  AND  DEVELOPMENT  CENTER, 
US  ARMY  AMCCOM ,  ATTN:  BENET  WEAPONS  LABORATORY,  SMCAR-LCB-TL , 
WATERVLIET,  NY  12189,  OF  ANY  ADDRESS  CHANGES. 


TECHNICAL  REPORT  EXTERNAL  DISTRIBUTION  LIST  (CONT’D) 


NO.  OF 

NO.  OF 

COPIES 

COPIES 

COMMANDER 

DIRECTOR 

US  ARMY  MATERIALS  &  MECHANICS 

US  NAVAL  RESEARCH  LAB 

RESEARCH  CENTER 

2 

ATTN:  DIR,  MECH  DIV 

1 

ATTN:  TECH  LIB  -  DRXMR-PL 

CODE  26-27,  (DOC  LIB) 

1 

WATERTOWN,  MA  01272 

WASHINGTON,  D.C.  20375 

COMMANDER 

COMMANDER 

US  ARMY  RESEARCH  OFFICE 

AIR  FORCE  ARMAMENT  LABORATORY 

ATTN:  CHIEF,.  IPO 

1 

ATTN:  AFATL/DLJ 

1 

P.O.  BOX  12211 

AFATL/DLJG 

1 

RESEARCH  TRIANGLE  PARK,  NC  27709 

EGLIN  AFB,  FL  32542 

COMMANDER 

METALS  &  CERAMICS  INFO  CTR 

US  ARMY  HARRY  DIAMOND  LAB 

BATTELLE  COLUMBUS  LAB 

1 

ATTN:  TECH  LIB 

1 

505  KING  AVENUE 

2800  POWDER  MILL  ROAD 

COLUMBUS,  OH  43201 

ADELPHIA,  MD  20783 

COMMANDER 

NAVAL  SURFACE  WEAPONS  CTR 
ATTN:  TECHNICAL  LIBRARY  1 

CODE  X212 

DAHLGREN,  VA  22448 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARMAMENT  RESEARCH  AND  DEVELOPMENT  CENTER, 
US  ARMY  AMCCOM,  ATTN:  BENET  WEAPONS  LABORATORY,  SMCAR-LCB-TL , 
WATERVLIET,  NY  L2I89,  OF  ANY  ADDRESS  CHANGES. 


